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Abstract. The long term perturbation of a Newtonian binary 
system by an incident gravitational wave is discussed in connection 
with the issue of gravitational ionization. The periodic orbits of 
the planar tidal equation are investigated and the conditions for 
their existence are presented. The possibility of ionization of a 
Keplerian orbit via gravitational radiation is discussed. 



1. Introduction 

A Newtonian two-body system cannot be completely isolated from 
all other masses in the universe as a consequence of the universality of 
the gravitational interaction. In fact, the attraction of the other masses 
would cause the binary system to move through approximately inertial 
spacetime. This center-of-mass motion should be distinguished from 
the relative motion, which is affected by the gradient of the disturbing 
forces. Consider, for instance, the equations of motion for an "isolated" 
two-body system in Newtonian mechanics 

d^Xi Gomim2 „ . r-^^/^ x 

™.^ + ^^(x.-xo = -m.v*(x,), (1) 

where Go is Newton's constant of gravitation and 

Gonip 



$(X) :=-E 



X — Xpl 
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represents the combined gravitational potential of all other masses nip 
at Xp in the universe. Here and throughout this work, the finite size of 
astronomical bodies is neglected. If, in the inertial space coordinates 
(X^,X^,X^), the binary system is so far away from the other masses 
that the relative distance between the masses comprising the binary 
is very small compared to the distance of the center of mass of the 
binary to the external masses, then, to first order in this small ratio, 
the equation of relative motion has the form 

d^r* kr 



^+^--'<.Mr', (2) 

where r = (r^, r^, r^) := Xi — X2, r is the length of r, and k = Go{mi + 
7712). Here, Kij, the tidal matrix, is given by 

9^* 



dX^dXi 



evaluated at the center of mass of the binary system. In Newtonian me- 
chanics the gravitational potential $ is a harmonic function; therefore, 
the symmetric tidal matrix is trace-free. 

It turns out that (0) holds approximately in general relativity as well, 
except that Kij would be represented by the "electric" components of 
the Riemannian curvature of the underlying spacetime projected onto 
a Fermi frame along the center-of-mass worldline [^, |I0|. That is, the 



equation of relative motion can be considered to be the Newton- Jacobi 
equation in the sense that once the internal Newtonian attraction is 
neglected, equation (g) reduces to the Jacobi equation in Fermi normal 
coordinates for the relative motion of two neighboring geodesies in the 
underlying spacetime manifold. Thus the spacetime coordinates in (Q) 
refer to a local Fermi system established along the path of the center 
of mass of the system. In our approximate treatment, we neglect rela- 
tivistic effects in the binary system. On the other hand, the external 
infiuences may now include gravitational radiation. It should be noted 
in this connection that classical celestial mechanics has been mainly 
concerned with the n-body problem; however, the "vacuum" between 
these bodies is expected to abound with gravitational radiation as well 
as with other radiation fields. It is therefore interesting to consider 
the interaction of gravitational waves with n-body systems, since it 
is estimated that half of all stars are members of binary or multiple 
systems. 

In this paper, attention is focused on a Newtonian binary system that 
undergoes perturbation due to an incident gravitational wave. Let the 
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spacetime metric due to the gravitational wave be given by 

where rj^y is the Minkowski metric, e is the strength of the perturbation, 
< e ^ 1, and Xtiu represents the gravitational radiation field. In the 
transverse traceless gauge, Xo^i = and Xij is a symmetric traceless 
matrix that satisfies the wave equation o'^Xij = and the transversality 
condition djXij = 0. It turns out that in this gauge, 

K.,{t) = -le^{t,X^^), (3) 

where (mi + m2)Xcm = miXi + 7712X2- It is possible to fix the posi- 
tion of the center of mass (e.g., Xcm = 0) in the approximation un- 
der consideration here, since Kij is considered only to first order in e. 
The perturbing field Xij be expressed as a Fourier sum of plane 
monochromatic waves with wavelengths much larger than the semima- 
jor axis of the binary system. Such waves could be generated by the 
motion of masses during the Hubble expansion, or could be primordial 
waves left over from the big bang era. It is therefore important to note 
that in our analysis there is no need to specify the initial conditions for 
the interaction of the waves with the binary; instead, we concentrate 
on the "steady-state" situation involving a dominant plane wave of fre- 
quency Q incident on the binary. Hence, the symmetric and traceless 
tidal matrix in equation (^ is given in Cartesian coordinates by 

Kii =eafi^cos^ 9 cos (fit), 
K12 =eP^l^ cos 6 cos (fit + p), 

= — eafi^ cos O sin cos (fit) , 
K22 = — eafi^ cos (fit), 
7^23 = - e/5fi^ sin 9 cos (fit + p), 

where a, (3, and p are constants, and O is the polar angle from the 
normal to the plane of the unperturbed orbit to the propagation vector 
of the incident radiation. However, for the sake of simplicity we will 
explicitly consider only the case of normal incidence, i.e. G = 0. In 
this case, the orbital plane will remain fixed; that is, when the waves 
are normally incident, the problem reduces to planar motion under the 
approximations considered here. 

Let the unperturbed Keplerian motion be confined to the (x, y)- 
plane. The transverse nature of the radiation field implies that the 
orbital plane will be unchanged under the perturbation of normally 
incident waves. Thus in (^, we have = x, = y, and = 0. The 
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nonzero elements of the tidal matrix for our single-frequency radiation 
are 

Kii =—K22 =eafi^ cos 

Ki2 = K21 =e/3n^ cos {nt + p), 

where ea and e/? represent the amplitudes of the two independent lin- 
ear polarization states of the low-frequency gravitational wave, and p 
represents the constant phase difference between them. 

The justification for replacing the actual problem with this rather 
simplified nonlinear model is that it becomes amenable to mathemat- 
ical analysis. It should also be remarked that — within the limitations 
discussed in this section — equation (|^) for the relative motion holds 
generally in the Fermi coordinate system established along the center- 
of-mass worldline. Thus, for this system, Kij = Roioj, where R^upa 
denotes the Riemannian curvature due to external sources projected 
onto the frame of the center of mass. In the Newtonian limit of general 
relativity, each Kij reduces to a second partial derivative of the exter- 
nal Newtonian potential $(X) evaluated along the path of the center of 
mass. On the other hand, for a weak gravitational wave, equation 
holds to first order in the amplitude of the gravitational potential. 
Thus, in general, the matrix (Kij) is a function of the proper time 
along the path of the center of mass. It is always possible to diagonal- 
ize this symmetric matrix; however, its dependence upon time implies 
that (0) must then be written in a rotating system of coordinates. 

In electrodynamics, the interaction of electromagnetic radiation with 
a two-body system constitutes a basic problem (e.g., the scattering and 
absorption of light by the Rutherford-Bohr atom). The gravitational 
analog of this problem in the classical regime would involve the scatter- 
ing and absorption of gravitational radiation by a Keplerian two-body 
system. The wavelength of light is much larger than the Bohr radius of 
the atom; therefore, the dominant interaction takes place via the elec- 
tric dipole moment of the atom since electromagnetism is a spin - 1 field. 
We expect by analogy that for gravitational radiation with a (reduced) 
wavelength that is much larger than the semimajor axis of the Keple- 
rian orbit, the dominant interaction would involve the mass quadrupole 
moment of the binary since gravitation is a spin - 2 field. This approx- 
imation corresponds precisely to dropping higher-order terms in the 
tidal equation (||). 

The reciprocity between emission and absorption of radiation should 
be noted. In the quadrupole approximation for the emission of gravita- 
tional radiation, the waves carry away energy and angular momentum 
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but not linear momentum. The same holds in the inverse process as 
well, except that in general the system can gain or lose energy. More- 
over, a Newtonian binary system emits gravitational radiation of fre- 
quency fl = mu, m = 1,2,3..., where u is the Keplerian frequency 
of the elliptical orbit. Similarly, resonant absorption of gravitational 
waves by an elliptical binary occurs at = muj, m = 1, 2, 3 ... , ac- 
cording to the linear perturbation analysis |]10 



A two-body system continuously emits gravitational radiation ac- 
cording to general relativity. Gravitational energy in the form of radi- 
ation is thus carried away from the system. Hence, the relative orbit 
evolves in such a way that the semimajor axis of the osculating ellipse 
monotonically shrinks. This phenomenon of inward spiraling of the 
members of the binary is consistent with the timing observations of 
the Hulse- Taylor binary pulsar [0 |T^. Direct observational evidence 



for gravitational radiation does not exist at present; however, efforts 
are under way to detect in the laboratory gravitational waves emitted 
by astrophysical sources with e ~ 10~^°. 

In this work, we will ignore the emission of gravitational radiation 
by the binary system, and concentrate our attention instead on the 
absorption process. The flow of energy between the incident radiation 
and the binary is not unidirectional, however. The self-gravitating 
system can absorb energy from the radiation field or deposit energy into 
the wave so as to induce an amplification of the radiation. These issues 
were first discussed in connection with the problem of ionization |10] in 
the context of linear perturbation analysis that in general breaks down 
over time. Here we employ the concepts introduced by Poincare |T^ 
for the treatment of nonlinear problems. These enable us to prove that 
periodic orbits exist in the perturbed system for which energy must 
steadily flow back and forth between the wave and the binary. It is 
important to emphasize that such periodic orbits occur near resonance 
conditions when certain definite phase relationships are satisfied. If the 
binary system monotonically absorbs energy from the wave, then the 
semimajor axis of the osculating ellipse will grow with time and the 
system eventually ionizes. We provide a qualitative picture for such a 
process in § ||. 

The gravitational quadrupole interaction may be illustrated by con- 
sidering the Hamiltonian for the relative motion 



n=^p'-'^ + ^K,,{t)Q,„ (4) 
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where the quadrupole moment per unit mass is defined by Qij — 
3r^r^ — r^5ij and in the most general case of a normally incident gravita- 
tional wave packet considered in this paper the matrix Kij is a traceless 
symmetric matrix of periodic functions. Thus Ku = —K22 = hi(t), 
K12 = K21 = h2{t), and there exist ri and T2 such that hi{t) = hi{t+Ti) 
and h2{t) — h2{t + T2). Here hi and /i2 represent the amplitudes of the 
two linearly independent polarization states of the perpendicularly in- 
cident gravitational wave. Now let 



denote the Newtonian energy, orbital angular momentum and the Runge- 
Lenz vector associated with the relative motion of the system (per unit 
reduced mass); then, these otherwise conserved quantities vary as a 
consequence of the coupling of the quadrupole moment of the system 
to the curvature of the background spacetime. Thus, in this quadrupole 
approximation, the Keplerian orbit exchanges energy and angular mo- 
mentum with the radiation field. At each instant, the relative mo- 
tion can be described by the orbital elements of the osculating ellipse. 
This osculating ellipse continuously makes transitions to other osculat- 
ing Keplerian orbits with different energies and angular momenta as 
a consequence of interaction with the external wave. It is interesting 
to describe the path of the system in the six-dimensional manifold of 
orbital elements; in fact, this paper is devoted to the description of 
periodic paths in this manifold. 

The interaction of gravitational radiation with matter may have 
played a significant role in the evolution of the universe. The treatment 
presented here is confined, however, to the interaction of an incident 
wave with a Newtonian binary system. In particular, we neglect all 
relativistic effects in the relative motion of the binary system. Let the 
system have a Keplerian frequency cu and semimajor axis a; then, by 
Kepler's third law, c<j^ = k/a^. Relativistic two-body effects may be 
neglected provided k <^ c^a, where c is the speed of light in vacuum. 
Moreover, the quadrupole approximation for the interaction of the sys- 
tem with the gravitational wave is valid if Jla <^ c. More generally, our 
approach is sound provided 



Furthermore, the requirement that the external wave be a small per- 
turbation of the system implies that fl/uj <C since the strength 
of the interaction is given by eO^/o;^, a quantity that must be much 
smaller than unity. Our objective is to determine the conditions under 
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which periodic Keplerian motions of the binary are continued to pe- 
riodic motions under the interaction. It is an important consequence 
of our methods, which are originally due to Poincare [Q, that the ex- 



istence of higher-order perturbing influences on the orbit, i.e., terms 
of order at least e^, can only affect the shape of the resulting periodic 
orbit but not its existence. In this first treatment of the nonlinear case, 
we consider only special cases of the other interesting questions — such 
as gravitational ionization — that are suggested by the electrodynamic 
analogy discussed in |]rU|. A treatment of the general problem on the 
basis of linear perturbation theory is contained in previous work |TD 



Superposition of linear perturbations due to the Fourier components 
of a pulse of gravitational radiation permits a general analysis in that 
case; however, the validity of the treatment is restricted in time as tem- 
poral evolution leads to a breakdown of the linear perturbation theory. 
Therefore, the intrinsic nonlinearity of the system must in general be 
taken into account for applications in celestial mechanics. In this re- 
gard, we mention that the equations of motion of a binary influenced 
by the gravitational attraction of a massive distant third body can also 
be treated using the methods developed here; this constitutes a special 
limiting case of the three-body problem and is discussed in Appendix 
B. 

We will be mainly concerned with the continuation of Keplerian or- 
bits under perturbation by a resonant gravitational wave. In general, 
we show that all but a finite set of such resonant orbits are not con- 
tinued to periodic orbits under the influence of the incident wave and 
that in general all elements of the finite exceptional set are, in fact, 
continued. 

The plan of this paper is the following. In § |^, we transform the 
perturbation problem to Delaunay elements and obtain explicit ex- 
pressions for the transformed perturbation in terms of Fourier series 
with coefficients that involve the Bessel functions. In § 0, we outline a 
continuation method based on the Lyapunov-Schmidt reduction that is 
adapted from P]. In § the results of the continuation theory are ap- 
plied to the perturbation of a binary influenced by a normally incident 
wave. In particular, we find bifurcation equations for the problem, that 
is, a system of equations whose simple zeros correspond to the continu- 
able periodic Keplerian orbits and we show that these equations indeed 
have simple zeros. In § ^, we consider the special case of circularly po- 
larized gravitational waves, a case that is not covered by the results of 
§ ^ In this case we show that there are periodic solutions. We also 
show for sufficiently weak perturbations of Keplerian ellipses that the 
semimajor axis of the osculating ellipse remains bounded for all time. 
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The final section, § ^, contains a brief discussion of some additional 
speculative results on the ionization problem. Some standard formulas 
are relegated to Appendix A, and in Appendix B we consider a spe- 
cial case of the three body problem: A binary influenced by a distant 
massive third body. 

2. Delaunay Elements and Fourier Series Expansion 

In terms of the canonical variables {pr,Pd,r,6), that is 

X = r cos 9, y = r sin 9, Pr = r, pe = r'^9, 

the Hamiltonian for our perturbation problem, with perturbation pa- 
rameter e, may be expressed in the following general form 

^' " ^ + er^[(p[t)cos29 + f[t)sm29), 

(5) 



n{pr,pe, r, 9) = ^{pI + ^) - - + er\(t){t) cos29 + ^(t) sin 20), 



where and ^^ are periodic functions with a common period. We will 
continue to use this general form in order to show how our theory can 
be applied. However, for the sinusoidal monochromatic gravitational 
wave model we will consider only the case where 

0(t) = cos(fit) , V'(^) = cos(fit + p) . (6) 

The unperturbed Hamiltonian 

Hipr,po,r,9) = ^ipl + ^)-^ (7) 

is called the Kepler Hamiltonian. We will consider only those motions 
corresponding to negative energy E = H{pr,p0,r,9). 

Following S. Sternberg fl^. Vol. 2, pp. 234-247], we define 

and we let a(l ± e) denote the roots of the quadratic polynomial 

2 2L2 G^L^ 

so that 

a = y , e = 1(L^ - (8) 

Here, a is the semimajor axis and e is the eccentricity of the Keplerian 
ellipse with < e < 1. However, we will only consider the case e > 0, 
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that is noncircular orbits. With this restriction in force, we define m, the 
eccentric anomaly and the true anomaly, implicitly by the formulas 

ail — e^) 

r = all — ecosMi, r = , (9) 

^ ' 1 + ecosw ^ ' 

and new variables ^ and g by 

^ = u — e sin m, g = 6 — v. 



As proved in [|T3|, the change of coordinates 

{Pr,pe,r,e) ^ {L,G,i,g) 

is canonical. Here i and g are "angle variables", defined modulo 27r, 
while L and G are "action variables" . The new coordinates (L, G, i, g) 
are called the Delaunay elements. 

In Delaunay variables, our Hamiltonian (^) is transformed to 

n = + e(C(L, G, i, g)<P{t) + G, 

where C (resp. iS) is the function obtained by expressing cos 29 (resp. 

sin 29) in terms of the Delaunay elements. Using the fact that the 
change of coordinates is canonical, the differential equations of motion 
are given by the Hamiltonian system 

dC (99 

L = -<g^{L, G, i, g)m + — (L, G, £, g)^), 

dC dS 
G = -<g^iL, G, i, g)4>{t) + — (L, G, i, g)m), 

dC (99 

^ ='^+<aL ^' ^' ^^^^^^ + dZ^^' ^' ^' 

(9C (99 

^ = 6(— (L, G, £, + — (L, G, £, (10) 

where u := k'^/L^ is the frequency of the elliptical Keplerian orbit. 

In order to analyze system (10), we must find computable expressions 
for the partial derivatives of C and S. This can be done in several ways; 
however, for our purposes, the most useful expressions are obtained 
from Fourier series expanded as functions of the angle variable i. The 
determination of these series is outlined in Appendix A, and the result 
can be expressed as 

5 °° 
C(L, G, i, g) =-a^e^ cos 2g + a? ^(^4,^ cos 2(7 cos z/£ — sin 2g sin z/£), 

5 °° 
5(L, G, £, g) =-a^e^ sin + ^^'^ ^5' "^^^ + -^'^ cos 2(7 sin z/£), 

2 -1 (11) 
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where 

A, =^(2z/e(l - e^)Jl{ve) - (2 - e^) J.(z/e)), 

B- = - -^^^-^ {^J>^) - ^(1 - e')J,{ve)). (12) 

3. Continuation Theory 

In order to analyze the continuation (persistence) of periodic solu- 
tions of the Kepler system to system (10), we use a method proposed 
in Here, we outline the main ideas; the reader is referred to |^ for 
the details. 

Consider a system of the form 

M = F(m) + e/i(M,t), (13) 

where m is a coordinate on a manifold M consisting of a cross product 
of Euclidean spaces and tori, h is 2tt/VL periodic in its second variable, 
and e is a small parameter. Let t i— > u{t, ^, e) denote the solution of (|13D 
with initial condition u(0,^,e) = ^, ^ G M. Also, we define the mth 
order Poincare map by 'P'"(^,e) = 'u(2m7r/f2, ^, e); it corresponds to a 
strobe that illuminates the orbit after m cycles of the perturbation. Of 
course, a fixed point of ^ ^— 'P™'(^,e) corresponds to a periodic orbit 
of (|13|). If m is the smallest such integer for which ^ is a fixed point, 
then ^ is the initial point of a subharmonic of order m. 

Suppose that there is a submanifold Z <Z M consisting entirely of 
fixed points of the unperturbed order m Poincare map, defined by 
P^iC) •= ^"^(^^O). Our continuation theory is a method, one among 
many, to decide if any of these fixed points survive after perturbation. 
More precisely, we say a point z E and therefore the unperturbed 
periodic orbit of (|T^) with initial point z, is continuable (or that it 
persists) if there is a continuous curve e ^ 7(e) in M such that 7(0) = z 
and 'P''"(7(e), e) = 7(e). Here, 7(e) G M is the initial point of a periodic 
solution of (p^). 

In order to apply the method of 0, namely Lyapunov-Schmidt re- 
duction to the Implicit Function Theorem, the fixed-point manifold 
(resonance manifold) Z must satisfy a nondegeneracy condition rel- 
ative to the unperturbed Poincare map. To specify this condition, 
consider z E Z and the derivative Dp^{z) viewed as a linear trans- 
formation of the tangent space T^M. The base point stays fixed be- 
cause is the identity on Z. Moreover, every vector in T^M that 
is tangent to the submanifold Z is fixed by Dp'^{z), or, as we will 
say, every such vector is in the kernel of the infinitesimal displacement 
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= Dp"^{z) — I. The manifold Z is called normally nondegener- 
ate if the kernel of the infinitesimal displacement is exactly the tangent 
space TzZ C M. Equivalently, Z is normally nondegenerate, if for each 
z & Z, the dimension of the kernel of the infinitesimal displacement at 
z is equal to the dimension of the manifold Z. 

Suppose Z has dimension A and that it is a normally nondegenerate 
submanifold of M. In this case the range of the infinitesimal displace- 
ment at each point in Z has codimension A. Thus, for z & Z, there 
is a vector space complement S{z), to the range of V^z). We let s{z) 
denote the projection of T^M to S{z). By choosing local coordinates, 
we note that both Z and S{z) may be identified with M^. 

Let z & Z and consider the curve in M given by e i-^ V^{z, e). This 
curve passes through z at e = 0. Its tangent vector at e = 0, which may 
be identified with the partial derivative V^{z, 0), is in T^M. We define 
the bifurcation function B to be the map, from Z to the complement 
S of the range of the infinitesimal displacement, given by 

B{z) = ~s{z)Vr{z,0). 

In local coordinates 8 -.R^ -^R'^. We will say z G Z is a simple zero 
of the bifurcation function provided B{z) = and the derivative DB{z) 
is invertible. 

A result in |Q is the following continuation theorem: 

Theorem 3.1. If Z is a normally nondegenerate fixed-point submani- 
fold of M for system (|13]) and if z & Z is a simple zero of the associated 
bifurcation function, then the unperturbed periodic orbit of (|13]) with 
initial point z is continuable. 

To use Theorem ^]T] as a practical tool, we must be able to compute 
V^{z, 0). Fortunately, this partial derivative can usually be computed. 
In fact, if n = fi(e), then 

0) = _-^r]'(0)M(2m7r/l](0), z, 0) + M,(2m7r/fi(0), z, 0). 

Thus, if t I— ^ W{t) is the solution of the second variational initial value 
problem 

W = DF{u{t, z, 0))W + h{u{t, z, 0), t), W{fd) = 0, 

then 

2^77 TT 

Vriz, 0) = -^n\0)F{u{2mn/QiO), z, 0)) + Wi2m7r/n{0)). 
In effect, W{t) = Ue{t, z, 0) with W{0) = because u{0, z, e) = z. 
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In our gravitational radiation model, it seems appropriate that the 
frequency of the gravitational wave is independent of the amplitude of 
the wave. Thus, we will assume below that fl does not depend on e. 
This simplifies the expression for V^{z, 0) by removing the "detuning" 



4. Continuation of Kepler Orbits 

To apply Theorem |3.1| to our perturbation problem (10), we must 
define a normally nondegenerate fixed-point manifold. For this, we 
consider the Kepler orbits that are in resonance with the periodic per- 
turbation. 

If there are fixed relatively prime positive integers m and n and a 
fixed value of u, the frequency of the Keplerian orbit, such that 

27T 271 

m— = n — , 

then the unperturbed solution of (10) starting at {L,G,£,g) is given 
by 

t {L,G,ujt + i,g), 

where t = t — to and to is an integration constant that denotes the 
starting instant of time. A detailed analysis shows that to can be set 
equal to zero here without loss of generality. Since £ is defined modulo 
2tt, this solution is periodic of period 27t/uj. Moreover, the mth order 
unperturbed Poincare map is defined by 

p"'mGJ,g) = mG,2iTm^ + e,g). 

If we define the three-dimensional manifold 

:= {{L,G,e,g) : moo = nfi}, 

and recall that £ and g are defined modulo 27r, then it follows immedi- 
ately that is fixed by p. To check that is normally nondegener- 
ate, we compute 



V{L,G,i,g) = Dp^{L,G,i,g)-I 



( \ 



-Q-nn/L 

/ 
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and we note that the infinitesimal displacement has a three-dimensional 
kernel that is spanned by the usual basis vectors 
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Moreover, the range of the infinitesimal displacement is complemented 
by the span of the vectors 
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To compute the bifurcation function associated with (10), we must 
compute the partial derivative V^IG, L, g,£,(}) on the manifold 
and then project the result into the complement of the range of the 
infinitesimal displacement. To do this, we simply solve the variational 
initial value problem 

dC dS 
4 = - -^{L, G, ut + g)4>{t) - — (L, G, ut + g)m, 

= - Q^{L, G, ujt + e, g)ct>{t) - — (L, G, ujt + £, g)m. 

1.2 an 

4 —JT^e + ^(^' ^' + ^' 9)<Pit) + ^{L, G, ujt + £, g)m. 

dC dS 

9e = -Q^{L, G, ut + e, g)(t>{t) + — (L, G, + £, y)V'(i) 

with zero initial values and then project the solution computed at 
t — m27r/Q into the complement of the range of the infinitesimal dis- 
placement. Prom this procedure, we obtain the following bifurcation 
function 



B(G,^,5) = (S^(G,^,^),S«(G,^,5),S^(G,^,^)), 



where 



dX 



B^{G,i,g):=-^, B''{G,i,g) := B^{G, i, g) := ^ 



dl 



dl 



de 



dg 



dG, 



(14) 



and 



X:= [C{L,G,ujt + e,g)(l){t) + S{L,G,ujt + i,g)i;{t) 



dt. 
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Using the resonance relation, we have 

1= / C{L,G,—nt + i,g)(t){t) + S{L,G,—nt + i, 

Jo m m 

and, after changing the variable io a = Vtt / m + 1 / we obtain 



9)m 



dt 



m 



n Je/n 



C{L, G, na, g)(p{m{a - i/n)/n) 



+ S{L, G, na, g)ipim{a - i/n)/n) 



da. 



Using the fact that the last integrand is periodic with period 2tt as a 
function of a and substituting (p and given in (P), we find 



X 



aC{L, G, na, g) cos(ma 



ml 



+ PS{L, G, na, g) cos{ma h p) 



n 



da. 



To compute X, we substitute the Fourier series (11) for C and S into 
the last expression for X and use trigonometric relations together with 
the fact that m and n are relatively prime, to conclude that X = 
unless n = 1. Of course, there may be continuable orbits for n > 1, 
but they are not detected by our first order method. In case n = 1, 
that is for the (m : n) = (m : 1) resonance, we find that 

X = -nma^Qi^a^Am cos m£ cos 2g — sin m£ sin 2g) 

+P{Am cos(m£ — p) sin 2(7 + BmSm{mi — p) cos 2g) 

This result can be rewritten as 
1 



X =-7rma^r2 
4 



{Am + -Bm)(tt — P sin p) cos(2(7 + mi) 

+ {Am — Bm) {a + P sin p) cos(2(7 — mi) + (3 {A^. + -Bm) cos p sm{2g + 

+ P{Am — Bm) COS p sm{2g — mi 

It is possible to express equation (15) in a more compact form, in the 
usual manner, by defining 

£ cos a := (3 cos p, £ sm. a := a + (3 sin p, 
jFcos r := /5cosp, jFsin r := a — /3 sinp, (16) 

so that 

£ = {a^ + + 2a(3smp)^^\ T= (a^ + /^^ - 2a/5sinp)^/^ 

and 

X= -nma^Q {Am+Bm)J^sm{2g+mi+T) + {Am—Bm)£sin{2g—mi+a) 
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The simple zeros of the bifurcation function are then the same as the 
simple zeros of 

T{Ara + S„) cos(2^ + + r) =0, 
S{Am — B„i) cos(2(7 — mi + a) =0, 

+ ^„ Lr sm(2(y( + mi + t) + — — — - \csm[2g — mi + a) =0. 



To obtain explicit formulas for the partial derivatives of the functions 
Ara and Bm with respect to G, we assume that G > so that 

dG ^ Le- 
In case G < 0, the partial derivative has the opposite sign and the 
subsequent analysis is similar. We use (32),(33),(|3^) (from Appendix 
A), and (12) to obtain 

^ = _ -±_[i2m\l - + 4)J„(me) - me(6 - e')J'^{me)), 

~d^^~ L^^( ~ ^""^^ " ^^)Jm{me) + e((2 - e^)(l - m^e') + m^)J'^{me)^ 



The simple zeros of the function B given by ([T^) correspond to the 
continuable periodic orbits. Equivalently, the continuable periodic or- 
bits correspond to the simple solutions of the system of equations given 
by (17). In order to find the simple zeros of (17), we will use the fol- 
lowing proposition: 

Proposition 4.1. If £T ^ (that is, \a\ ^ or \a\ = \(3\ but 
I sinpl < 1) and if the system of equations (17) has a solution, then 



(19) 



IS zero. 



Proof. If (17) has a solution and equation ([T9|) does not vanish, then, 
since the first factor of equation (|19]) is not zero, we have 

cos(25( + mi + T) = 0, cos(25( - mi + a) = 0. 

This implies that 

sm{2g + mi + T) = ±1, sm{2g - mi + a) = ±1, 

and, since the third equation of (17) is zero, that 



\dG dG J V dG dG 
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in contradiction to the fact that (p!9D is not zero. □ 



Proposition ^?T] reduces the search for solutions to several cases. Just 
note that as soon as one of the factors of ([T9|) vanishes, the value of e 
and hence the values of 

„ . „ dAm dBm 

4 _ R 4 _|_ R 

oG oG 

are fixed. Thus, equations (17) reduce to solvable trigonometric equa- 
tions. It is important to note that Proposition |4.1| does not cover 



the interesting case of circular polarization, which is therefore deferred 
to §|. 

To study the zeros of ([19|) we will use the following proposition. 

Proposition 4.2. For the (1 : 1) resonance, the functions Ai + Bi and 
Ai — Bi appearing in equations (17) and viewed as functions of e are 
both negative on the interval < e < 1. The functions dAi/dG and 
dBi/dG, viewed as functions of e, each have exactly one simple zero on 
the interval and their zeros are distinct. For the (m : 1) resonance with 
m > 1, the range of the function {dAm/dG)/{dBm/dG), viewed as a 
function of e on the interval < e < 1, contains the interval [—1, 1]. 

Proof. We will outline the proof of the proposition. Some of the com- 
putations were checked using a computer algebra system. 

Consider the case m = 1. We will use the following elementary 
lemma Lemma 3.5] to show that the function / defined by 

/(e) ■.=j{A,-B,) 



=2e((l - e^) + VT^)J[{e) - ((2 - e^) + 2(1 - e^)^/^) J,(e) 

(20) 

is negative on the interval Jq := {e : < e < 1}. 

Lemma 4.3. Suppose f is a smooth function defined on an interval 
[a, b) with the additional property that there is a number e > such 
that f{x)f'{x) > for a < x < a + e. If there are smooth functions p, 
q, and r defined on {a,b) such that p{x)r{x) > and 

p{x)f"{x)=q{x)f\x)+rix)f{x) (21) 



on the interval {a,b), then f is strictly monotone on [a,b). In particu- 
lar, f{x) has the same sign on (a, b) that it does on (a, a + e). 
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The function / defined by (20) satisfies a differential equation of tlie 
form (pT]) with x = w := — and 



p(e) 
g(e) 



re 



=e^w^(5«; + 2), 
=e^(15«;^ + W + 2), 

=6^(5^^ - - 16w^ + 26w^ + 2Aw + 4). 



To test the sign of p(e)r(e), we change the variable to w and note that 
< < 1. Let 



p*{w) := p{Vl — w"^), r*{w) := r(-\/l — 



2^ 



Clearly, is positive on < w < 1. The second factor of r* is easily 

shown to be positive on the same interval. For example, the second 

factor is positive at w = and has no roots in the interval. The fact 

that there are no roots can be checked by computing a Sturm sequence 

(cf. 0). This proves p(e)r(e) > for e G Jq. 

To complete the proof of this case, it remains only to show that 

there is some e > such that /(e) < and f'{e) < on the interval 

< e < e. To do this, we simply note that the Taylor series of / and 

/' at e = are given by 

7 35 
/(e) = + O(e^), /'(e) = -|e^ + 0(e«). 

The fact that Ai + Bi is negative on the interval Jq can be proved 
in a similar manner. 

We must show that dAi/dG has exactly one simple root on the 
interval Jq. To do this it suffices to prove that the function given by 

e H-^ (6 - e^) J2(e) + e{2e^ - 3) Ji(e) 

has exactly one simple zero on Jq. Equivalently, using the fact that 
J2(e) does not vanish on Iq, it suffices to show that the function 

has the same property. This fact follows from the expression (18) for 
dAi/dG, the recurrence formula 

jyJ,y{x) - xJl{x) = xJy+i{x), (22) 

that is a simple consequence of equation (32) of Appendix A, and the 
connection between e and G. 

The fact that /o has exactly one simple zero is a consequence of the 
following three propositions: (a) The function /i(e) := (6— e^)/(2e^— 3) 
is monotone decreasing on Iq. (6) The function /2(e) := eJi(e)/J2(e) 
is monotone decreasing on Iq. (c) The function /o has a zero in Iq. 
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Statement (a) is immediate: f[ is negative on the interval. Statement 
(6) follows from the product representation of the Bessel function J^, 
given by [0 

where j,y^s denotes the sth zero of and the fact that the zeros are 
interlaced as follows: 

ju,l < ju+1,1 < iu,2 < ju+1,2 < ju,3 < ■ ■ ■ 

Statement (c) follows from two facts: lime^o+ /o(e) = 2 and /o(l) < 0. 
The first fact is immediate from the above product representation or 
from the Taylor series for the Bessel functions at e = 0. To obtain the 
inequahty, we use the product representation of the Bessel functions to 
deduce 

OO -1 OO -1 

/o(i) = -5+4n(i--^)/n(i--^)- 

s=l s=l J2,s 

Since j2,s > ji,s > 1 for all s, the quotient of the two products is less 
than unity, as required. 

The fact that dBi/dG has a unique simple zero is proved using a 
similar analysis. Just as before, it suffices to show that the following 
function has exactly one simple zero on Jq: 

/(e) :=-e3:^ + e^-3e^ + 3. 

Both terms are monotone decreasing on Jq with /(O) = 3 and /(I) < 0. 

We claim the zeros of OAi/dG and dBi/dG do not occur at the same 
point. After division by nonzero factors and the substitution m = 1 
in (18), it suffices to show that the following functions do not have a 
common zero on Jq: 

S(e) : = (2(1 - e^f + 4) Ji(e) - e(6 - e^)J[{e), 

P(e) := - 3(1 - e^) Ji(e) + e((2 - e^){l - e^) + l)J[{e). 

If the functions do have a common zero, then the function 

/(e) := S(e) - e^(e) 

has at least one zero on Iq. To prove the claim, we show that / is 
negative on Iq. 

Using the recurrence formulas for Bessel's functions, we find that 
/(e) = (-e''+3e^+e^-3e2-6e)Jo(e) + (e^+2e^-6e^-5e2+6e+12)Ji(e). 
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To prove this function is negative on Jo we will apply Lemma [4.3|. We 



find that / satisfies a differential equation of the specified form with 



p{e) 
g(e) 
f(e) 



=e2(e^° - 7e^ - 4e^ + 31e^ + 14e^ - 55e^ - 456^ + 846^ + 9e - 36), 
=e(lle^° - 63e® - 32e^ + 217e^ + 84e^ - 275e^ - ISOe^ + 2526^ + 18e 
= - _ 27eio + 126^ + 95e* + lle^ - 235e^ 
+ 273e^ - 19e^ - 4686^ + 3516^ + 90e - 108. 



Using Sturm sequences, it can be proved that p(e)f(e) > for e G Iq. 
Moreover, we find that 

^3 ^ ^ 97 

/(e) = -^e' + O(e^), /(e)/'(e) = ^e^ + O(e^). 

This completes the proof of the claim. 

In case m > 1, it suffices to consider the range of the function 
given by e (— > {dAm/de)/{dBm/de). A computation shows that the 
Taylor series of both the numerator and the denominator of is 
given by — 5e + O(e^) in case m = 2 and, in case m > 2, by 

8m'"(m-l)(m-2) ^^_3 ^ 

It follows that lime^o+ -fm(e) = 1. 

We claim that e t— > dBm/de has at least one zero on Iq. If not, then 
Bm is a monotone function of e. A computation shows that Bm has a 
removable singularity at e = and that 

If m > 2, then 5^ is increasing for < e << 1. But, from the 
definition of Bm, we have lime^i- Bm{e) = 0, in contradiction. For 
m = 2, we find that 

52(e) = 1-^6^ + 0(6^) 

and B2 decreases for < e << 1. But, B2 is negative for << e < 1. 
To see this just note that near e = 1, the sign of B2 is determined by 
— J2(2e). By standard properties of the Bessel functions (cf. [|l|), J2{x) 
is positive on the interval (0, j2 1), where Jl^ denotes the sth zero of J^. 
Since z/ < Jl^, we have that — J2(2e) < for << e < 1. Again, since 
limg^i- B2{e) = 0, we have a contradiction. This proves the claim. 

Suppose for the moment that dBm/de > dAm/de on the interval 
< e < 1 and consider the first zero oi e ^ dBm/de. It follows 
that dAm/ de{e^) < while dBm/de > on the interval < e < e^,. 
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Thus, we have lim^^^- Fm{e) = — oo and the range of contains the 
interval {—oo, 1], as required. 

To complete the proof, it suffices to show that the function given 

by 

e ^— s> irn^e'idBra/ de — dAm/de^ 

for m > 1 is positive on the interval < e < 1. This fact follows from 
Lemma |4.3| . The Taylor series of Gm at e = is given by 

2™"m\{m + 1) 

Thus, it follows that Gmie)G'^{e) > for < e << 1. We also find that 
there are functions Pm(e), 5'm(e) and r^ie) such that Pm(e)rm(e) > 
and 

Pm{e)G'^{e) = g^(e)G^(e) +r„(e)G„(e) 
on the interval < e < 1. In fact, 



Pm ■.=e^w^(5w^m^ + {-8w'^ + 30w^)m^ + (4w^ + 24w)m + 8 

+ {12w^ + 24w^ + 48w)m + 24), 
Tm ■.=ew'^m(5w^^m^ + {-18w^ + 60w^)m^ + {-12w^ - 52w^ + 180w 
+ (30w^ - 76w^ + I66w^)m - I2w^ + 32w^ - 12^2 + 24 



(To verify that Gm satisfies the second order differential equation with 
these coefficients, we compute the derivatives of Gm and then convert 
all the expressions to the variable w.) Finally, to show Pm(e)rm(e) > 0, 
it suffices to show that the inequality holds for < w < 1. To do this, 
view Pm and as polynomials in m and note that all their coefficients 
are positive functions of w. □ 

4.1. The (1 : 1) Resonance. The fundamental physical result of this 
section is the following proposition: Among the periodic Keplerian 
orbits in (1 : 1) resonance with an incident gravitational wave, there 
are generally a (nonzero) finite number of continuable periodic motions. 
In fact, the frequency of the gravitational wave fixes the semimajor axis, 
a, while the amplitudes and the phase shift, a, /3, p, of the wave fix 
the eccentricity of the unperturbed Keplerian orbits that are excited 
by the perturbation. The inclination of the major axis and the angular 
position on the ellipse that complete the set of initial conditions for a 
continuable orbit on the excited ellipse are given by formulas presented 
below. However, two facts complicate the mathematical analysis: there 
are exceptional choices of the wave amplitudes a and /? such that none 
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of the periodic orbits in (1 : 1) resonance with the incident gravitational 
wave are continuable and there are zeros of the bifurcation function that 
are not simple. 

The precise mathematical result that we will prove requires a gener- 
icity assumption. For this we will say that a property of the zero set 
of (17) is generic relative to the parameter vector p) G M^, if it 

holds for an open and dense subset of M^. 

Proposition 4.4. If m = 1, then, generically relative to the parame- 
ters {a,l3,p), the zero set of system (17) is a nonempty finite set con- 
sisting entirely of simple zeros. Ifm = 1, a^+(3'^ ^ 0, and ajS sin p = 0, 
then system (17) has a nonzero finite number of zeros which are all 
simple. 

Proof. The first generic assumption is SJ-" ^ 0, the second generic 
assumption is that o? ^fp' ^ 0. (Of course, if + = 0, then there is 



no perturbation of the Keplerian orbits.) According to Proposition [4.2| , 



we have — 7^ for < e < 1. Thus, if system (17) has a solution. 



then, according to Proposition [4.1| , we must have 



Define 



dG dG> ^dC dG 



and note that, after a simple algebraic manipulation and after taking 
into account the obvious fact that the partial derivatives with respect 
to G can be replaced with no loss of generality by partial derivatives 
with respect to the eccentricity e, the last condition is equivalent to 
the requirement that 

^ de de ■'^ de de ' 

Also, taking into account the fact that £ >0 and > 0, our assump- 
tion that + 7^ implies < |k| < 1. 

To find the solutions of system (17), suppose for the moment that 
the equation 

(23) 

de " de ^^^^ 

has a solution. For this value of e the third equation of system (17) 
vanishes provided sm{2g + £ + r) and sm{2g — £ -\- a) are both equal 
to one or both equal to minus one. In either case, cos(2(yf + £ + r) and 
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cos(25f — ^ + a) both vanish. Thus, for all integers M. and M such that 

2c/ + £ + r = ^ + 27rA^, 2g - i + a = ^ + 2'kM , 
or such that 

2g + l + T = -'^ + 2txM, 2g-i + a = -^ + 2nX, 

the fixed value of e together with the (nonzero) finite number of si- 
multaneous solutions of these last equations with the property that 
< f?,^ < 27r give a set of solutions of system (17). A similar result is 
valid in case 

dBi dAi 

To determine the simplicity of these solutions, we must compute the 
Jacobian of system (17) with respect to the variables {G,i,g) at the 
given solution. This Jacobian is easily computed by expanding along 
the first column of the Jacobian matrix. Up to a nonzero constant 
multiple, we find the value of the Jacobian to be 

(•4;-BS((^+|^)^-(29+f+r)+(|^-^)fsm(2s-m)). 

In particular, using the fact that e is a monotone function of G, it 
follows that the solution {G,i,g) of system (17) is simple provided the 
corresponding solution e of equation (p3D, respectively (pl[), is simple. 

To finish the proof, we must determine the existence and simplicity 
of the solutions of equations ( P3| ) and (pij). 



If either a = 0, /? = 0, or sinp = 0, then k = and both equa- 
tions ( P^D and ( P^ ) have unique simple solutions by Proposition ^ 
This proves the second statement of the proposition. 



Since the left-hand sides of equations (123|) and (124| ) are both analytic 
functions of e, there are at most a finite number of solutions on the 
interval < e < 1. Moreover, since the map (a,/3, p) t— > /t(a,/3,p) is 
regular on an open and dense subset of M^, if some of the solutions of 
one of the equations are not simple, then there is an arbitrarily small 
perturbation of the triplet (a,/3,p) such that the perturbed equations 
have a finite number of simple zeros. 

The existence part of the first statement of the proposition is a con- 
sequence of the following facts. If k 7^ 1 (a generic assumption), then 
equation (^) has at least one zero. To see this, we note that the func- 
tion e ^— s> dBi/de has value —3 at e = and has limit 00 as e — ^ 1~ 
while the function e ^— dAi/de has value —3 at e = and has a finite 
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value at e = 1. As long as k 7^ 1, then the function 

dBi dAi 
de de 

has a negative value at e = and has limit 00 as e — > 1^. □ 

The proposition does not give a complete description of the zero 
set of system (17). However, since this description is reduced to an 
investigation of equations (|23| ) and (^4[) that are algebraic combinations 



of Bessel's functions, numerical approximations suggest the following 
scenario. If k 7^ 1, then the function 

dB^ Ml 

has exactly one simple zero on the interval 0<e<l. IfK=l, 
then (|25|) vanishes at e = and increases monotonically to 00 as e — 
1~. If K < 0, then the function 

OA, dB, 

has exactly one simple zero on the interval < e < 1. There is a 
number ^ 0.036 such that if < k < k*, then p6|) has exactly 
two simple zeros, while if k, > k*, then (^) has no zeros. If k = k^, 
then (p6D has exactly one zero which is not simple. 



Remark 4.5. In case (3 = 0, that is the wave is plane polarized in 
a very special direction, it appears that the zeros of ([T9|) are all near 
6=1. The smallest occurs for the (m : 1) = (1 : 1) resonance in case 
dAi/dG = 0. Its root is larger than e = 0.68. The root is larger for 
the higher order resonances. This suggests that only some "comets" 
could remain periodic after perturbation by a gravitational wave with 
this particular polarization. 

4.2. The (m : 1) Resonance, m > 1. For the (m : 1) resonance we 
will prove that there are perturbed periodic solutions for the generic 
a, (3 and p. This is the content of the next proposition. 

Proposition 4.6. If m > 1, then generically, relative to the parame- 
ters a, (3 and p, system (17) has at least one simple zero. 

Proof. It suffices to consider a, [3 and p such that EJ-' 7^ 0. Let g and 
I denote a solution of the equations 



251 + m£ + r = 7r/2, 2g - mi + a = tt/2, 
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and note that with this choice of g and £, system (17) has a solution 
provided G , equivalently e, is chosen such that 



dG dG ^ ^ dG dG 
Equivalently, as in Proposition ( |4.4| ), there is a solution provided 



dG dG ~ ^^^^ 

Since \k\ < 1, an application of Proposition ( [4.2[ ) shows that (^) has at 
least one solution. Moreover, since A"^ — is not the zero function, 
it has only a finite number of zeros for < e < 1. Also, its zeros 
do not depend on the choice of the parameters a, /3 and p. Thus, if 
necessary, after a perturbation of the parameters we can be sure that 
our solution of ( jn\ ) is not a zero of Af^ — B"^ and that it is a simple zero 
of the left hand side of (|27|) . As in Proposition (|4.4| ), it follows that the 



corresponding choice of {G,i,g) is a simple zero of system (17). □ 



5. Circularly Polarized Waves 

In this section we consider the equations of motion (^ for the case 
of a circularly polarized incident wave. This corresponds to the special 
case where, in the components of the tidal matrix K in (|^), we take 
a = (3 and p = ±7r/2. The minus sign corresponds to a right circu- 
larly polarized wave, while the plus sign corresponds to a left circularly 
polarized wave. We note that this is the main case excluded from 
the analysis of the previous section. There, the bifurcation function 
does not have simple zeros for the bifurcation problem corresponding 
to circular polarization. 

The equations of motion for the right circularly polarized wave have 
the form 

——r H 5" = — eafi^ (x cos ^It + y sin fit) , 

—TT H — T =eafi^(— a; sin ilt + y cos fit). 
dt^ r-^ 

This system can be treated in a similar manner as the analogous sys- 
tem that arises in Hill's lunar theory (cf. [12| |13|). The key idea is 



to view the system in a new Cartesian coordinate system that rotates 
relative to the inertial system with half the frequency of the gravita- 
tional wave. This factor of 1/2 is due to the fact that the wave has 
helicity +2. In these rotating coordinates, that we again call x and y, 
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Figure 1. Orbits of a Poincare map for the Hamil- 



tonian system with Hamiltonian (^). The parame- 
ters are a = l, Q = l, k = l, e = (left hand 
panel), e = .026 (right hand panel), and the energy is 
H{pr, Pe, r, 6) = H{0, 1, 1, 0). For post script figures con- 
tact carmenchicone.cs.missouri.edu 

we obtain the following system 

d'^x dy 1 2 , ^a; 

d'^y ^dx ,1 ky , ^ 

_| + S,--(j + ™)tfy+^ = 0. (28) 

We note that the replacements t —t and Q —Q leave the system 
invariant. Thus, it suffices to consider the equations of motion for 
either state of circular polarization. 

A remarkable fact, also utilized by Hill, is that system (28) is equiv- 
alent to a Hamiltonian system with Hamiltonian 

H = liX' + Y') + ^{yX-xY)-^+ '-aQ\x' - y'), 

where X := x — fly/2, and Y := y + Qx/2. This Hamiltonian is given 
in polar coordinates by 



H = lipl + $)-^-^Pe+ '-anV^cos29, (29) 



where Pr = {xX + yY) /r and pg = xY — yX, and in Delaunay elements 
by 

H = - |g + '-an'cm G, e, g). 

The Delaunay form of the Hamiltonian is expressed in action-angle 
variables and is in the correct form to apply the Kolmogorov-Arnold- 
Moser (KAM) theory (see for example |^ |^). Here, the Hamiltonian 



is degenerate. But, as pointed out in Sternberg |jT3|, Vol. 2, p. 257], the 
Hamiltonian system with Hamiltonian H' has the same trajectories as 
the original system, only the speed along trajectories is changed. More- 
over, the unperturbed part of H' satisfies the nondegeneracy assump- 
tion for the KAM theorem — its Hessian, with respect to the actions, 
has a nonzero determinant. Thus, the perturbed trajectory remains 
bounded in time, being trapped between two-dimensional invariant tori 
in the three-dimensional energy surfaces of our two- degree of freedom 
Hamiltonian. Thus, sufficiently weak circularly polarized gravitational 
waves do not "ionize" the Keplerian ellipses; that is, the osculating 
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semimajor axes do not become unbounded. This is illustrated in Fig. |l|, 
where "phase portraits" for a typical Poincare map for the Hamilton- 
ian system corresponding to (pQ]) is depicted. After an energy Hq is 
fixed, each orbit on the graph is produced by first choosing an initial 
point {pr, r) in the depicted plane and then by marking the position of 
the {pr, r) coordinates of the Hamiltonian orbit, with initial condition 
{pr, r), 6' = 0, and with po the implicit solution of H{pr,pe, r, 0) = Hq, 
at each time when 6{t) is a multiple of 27r and 6{t)6{0) > 0. In the 
actual computation, 6 is reset to zero each time a crossing is marked. 
The figure contrasts the foliation by invariant tori for the unperturbed 
system, where there appears to be an incidental resonance of order two 
and one of order three, with the existence of several invariant tori, as 
well as a strongly stochastic layer, for the perturbed Poincare map. 

We mention that the existence of periodic solutions of the equations 
of motion in the rotating coordinate system (these correspond to pe- 
riodic or quasiperiodic motions of the original system) can be proved 
along the lines of Poincare's periodic solutions of the first and the sec- 
ond kind for the restricted three-body problem. 

The continuation theory for the periodic solutions of the first kind 
does not depend on the perturbation terms, only on the Floquet mul- 
tipliers of the "circular" periodic orbits of the Kepler problem in the 
rotating coordinate system (see [|lT| [1^). The unperturbed orbits that 
continue are given by pe = C , Pr = 0, r = C'^/k for a fixed constant C. 

The continuation theory for the elliptical orbits, periodic orbits of the 
second kind, can be completed along the lines attributed to Poincare 
and subsequent authors as outlined in [^, Vol 2, p. 274]. However, 
these can also be found using the continuation theory of § ^. In the 
following brief outline of the procedure, we will use the Delaunay for- 
mulation of the equations of motion. 

After isoenergetic reduction, by implicitly solving for the angular 
momentum G in the perturbed Hamiltonian as a power series in e, and 
reformulation of the reduced system as a system of differential equa- 
tions with the timelike independent variable g, one obtains a system of 
two differential equations that are vr-periodic in g. The continuation 
theory of § ^ can be applied to this reduced system. 

Here, the Poincare section is given by the submanifold defined by 
g = 0, and the return map is an iterate of the strobe after each g in- 
terval of length vr. An {m : n) = {m : 1) resonant unperturbed orbit 
corresponds to an invariant one dimensional torus in the Poincare sec- 
tion. All such tori are normally nondegenerate due to the fact that the 
periods of the unperturbed orbits, in the reduced unperturbed system. 
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change monotonically with L (this is equivalent to the twist condi- 
tion for the Poincare map). The corresponding bifurcation function 
maps the angular variable ^ along the unperturbed orbit to the av- 
erage of the first order part of the reduced differential equation for 
the action L over the unperturbed resonant orbit with initial value 
^. In fact, the function is given (up to a nonzero constant multi- 
ple) by £ ^— s> (y4m(e) — -Bm(e))sinm£ where and Bm are defined 
in (12). This function has simple zeros (for almost all values of the 
eccentricity). Hence, the unperturbed resonant orbits have continua- 
tions. In particular, our method produces a periodic orbit of the form 
g t— > (L{g,e),i{g,e)) for the reduced system with independent variable 
9- 

Using the fact that G is implicitly given as a function of the form 

G:=GiLig,e)Jig,e),g,e), 

we see that G is also periodic in g. Finally, to obtain g as a function 
of time, we use the first order differential equation 

This last equation, at least for sufficiently small e, has solutions g{t) 
such that, for some period T(e) > 0, its solution satisfies g(t + T{e)) = 
g{t) — 2ti. Since g is an angular variable, the corresponding function 

t^{L{g{t),e),G{g{t),e),^{g{t),e),g{t)) 

produces a periodic solution of the original perturbation problem in the 
rotating coordinate system. These solutions are analogous to Poincare's 
periodic solutions of the second kind. 

6. Speculations, Conjectures and Numerics 

Will a Keplerian binary perturbed by an incident gravitational wave 
ionize? To make this question precise, we consider the unperturbed 
system to be a Keplerian ellipse, that is, the eccentricity e of the un- 
perturbed orbit satisfies < e < 1; equivalently, the energy of the 
unperturbed system defined by Hamiltonian (0) is negative. The cor- 
responding perturbed orbit (the Hamiltonian trajectory given by (|^)) 
generally does not lie on an ellipse. However, we define its osculating 
conic section at the instant the perturbed motion reaches the state 
{pr,pg,r,6) to be the conic obtained as the Keplerian orbit with this 
initial data, that is the Keplerian motion that would be obtained if 
the perturbation were "turned off" at this instant of time. To ionize, 
the flow of energy between the binary and the wave must turn uni- 
directional in a time averaged sense in the course of the perturbation 
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such that the binary would steadily absorb energy; in time, the binary 
system would be permanently disrupted and the two bodies would even- 
tually be infinitely far apart from each other. On the other hand, the 
basic equation of motion (H) breaks down once the semimajor axis of 
the osculating ellipse becomes comparable to the (reduced) wavelength 
of the incident gravitational wave. To study the route to ionization, we 
therefore introduce the notion of dissociation. We say the Keplerian 
ellipse determined by the initial data {pr,p0,r,6) at time to dissociates 
under the infiuence of the perturbation if at some later time the os- 
culating conic along the perturbed orbit is a hyperbola. Equivalently, 
if one wishes to remove the geometric language of this definition, the 
requirement for dissociation may be recast as follows: the Keplerian 
energy H {pr(t) , pg(t) , r (t) , 0(t)), where H is given by (^, defined along 
the perturbed orbit becomes positive in the course of time. 

The ionization question probably does not have a simple answer. 
However, two facts are clear. If the strength of the perturbation is suf- 
ficiently small, there are Keplerian binaries that do not ionize. Indepen- 
dent of the strength of the perturbation, there are Keplerian binaries 
that do dissociate. The first fact is proved in this paper: some of the 
resonant Keplerian orbits continue to periodic orbits of the perturbed 
system. We also recall that in the case of sufficiently weak circularly 
polarized incident gravitational waves, as discussed in § ^, none of the 
Keplerian orbits ionize. On the other hand, to see that dissociation is 
possible and to speculate on the fate of all orbits, we must review the 
geometry of our problem. 

Recall that Hamiltonian (^) defines a 2i-degree of freedom Hamil- 
tonian system. This system is equivalent to the three-degree of freedom 
Hamiltonian system given by 

n*iJ,Pr,Pe,s,r,e) = J + hpl + ^) - - + er2(0(s)cos2^^ + V^(s)sin2^), 

^ r r ^3Q^ 

where J is a "fictitious" action variable conjugate to the "time" s. Note 
here that the phase space of (|30| ) is six dimensional and that the five 
dimensional submanifold V given by 

i/(p„p,,^,r) = i(p^ + 4)-- = 

separates the phase space. This submanifold corresponds to the par- 
abolic Kepler orbits while a Keplerian binary corresponds to an ini- 
tial point in the region of the six dimensional phase space given by 
HipnVei f) < 0. Dissociation occurs provided the perturbed orbit of 
the Keplerian binary eventually crosses the manifold V. 
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To determine that some perturbed orbits do in fact cross the manifold 
P, in both directions, we simply compute the derivative of H along the 
perturbed orbit to obtain 

H = -2e(rp,.(0(s) cos 2^ + ^Ij{s) sin 26) + pe{^{s) cos 20 - (f){s) sin 20)). 

This derivative may be interpreted as a measure of the cosine of the 
angle between the perturbed Hamiltonian velocity field and the normal 
to the sub manifold V. It is apparent that there are open sets on V 
where H > 0. Thus, there are open subsets (obtained by reversing the 
flow on the boundary set) of the region H < such that every point of 
the subset corresponds to a Keplerian ellipse that dissociates. However, 
we emphasize the fact that time, represented by s in our three-degree 
of freedom system, is one of the variables under consideration when 
these open sets are determined. Thus, the initial data for a Keplerian 
ellipse that dissociates include an initial time t = to- We do not know, 
from this analysis, how far back in time the reversed orbits remain in 
the region where H < 0. 

A more sophisticated analysis might be based on the diffusion prop- 
erties of the orbits of H*. The geometric picture that is believed to 
hold for the dynamics of a nearly integrable Hamiltonian system with 
at least three degrees of freedom is easy to describe informally: there is 
a dense set of invariant tori coexisting with a dense set of orbits some 
of which are dense in their respective energy surfaces (see 0] |^). 

For the Hamiltonian ( pOD , it is easy to see that each five dimensional 
energy surface intersects V. In fact, each energy surface intersects the 
subsets of V defined by if > and H < 0. Thus, in the situation of the 
conjectured dynamics, for a sufficiently small perturbation strength, 
some of the orbits not on invariant tori ionize while a large subset of 
the orbits on invariant tori do not. Of course, some of the invariant 
tori of the perturbed system might intersect V; under our definition, 
the corresponding orbits will dissociate even though these same orbits 
will repeatedly return to the region where their osculating conies are 
ellipses. 

We note that the usual theory that is used to prove the existence 
of invariant tori for nearly integrable Hamiltonian systems, namely 
the KAM theory, is not directly applicable to the Hamiltonian given 
by (0) because the unperturbed system does not meet the required 
nondegeneracy conditions. In fact, this Hamiltonian is degenerate and 
isoenergetically degenerate (cf. 0, p. 408]). In particular, these facts 
are evident from the Delaunay action-angle coordinates for dS^) , where 
we see that the three-dimensional unperturbed invariant tori given by 
fixing J, L, and G do not even have dense trajectories because g = 0. 
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To obtain "nondegenerate" tori, one must consider the two dimensional 
tori given by fixing J, L, G, and g while leaving ^ and s free. Some of 
these tori may survive perturbation and, given their dimensions, it is 
possible that some of the perturbed tori are "whiskered" : they have sta- 
ble and unstable manifolds (each with dimension at most three). The 
existence of these invariant manifolds — together with the stable and 
unstable manifolds associated with periodic solutions (one dimensional 
invariant tori) and the intersections among them — is likely responsible 
for the diffusion of some of the orbits not in the union of these invariant 
sets and the invariant tori. 

Figure 2. Projections into {pr,r) plane of part of 
one orbit, approximately 5000 iterates in each panel, 
of the Poincare map for the Hamiltonian system with 
Hamiltonian (^. The parameters are a = f3 = 2, 
k = 1, p = 0, and e in the panels from left to right 
is 0.0, 0.001, 0.002 and 0.0025. The initial values are 
{Pr,Pe,T,9) = (.5,1,1,0). In this case, fl is chosen 
[fl ?a 3.897) so that the unperturbed Keplerian ellipse 
has frequency approximately l/6th the frequency of the 
incident gravitational wave. The region bounded by the 
branches of the curve given by rp^ = 2k shown in the 
panels corresponds to elliptical motion. To obtain the 
.ps files for this figure or hard copies of the figure, 
cont act carmenchicone . cs . missouri . edu 

We end this section with a short description of some of the numerical 
experiments performed on the Hamiltonian system (^) given by 

n{pr,pe,r, 9) = -{pI + ^)-- + er\<P{t) cos 20 + ^(t) sin 20), 

where and ip are given by equation (||). The results of a typical 
experiment that suggests the possibility of dissociation for an elliptical 
Keplerian orbit with eccentricity e = 0.5 are depicted in Fig. |^. To 
obtain the figure, the above 2|-degree of freedom Hamiltonian system 
is integrated numerically, and the values of the solution are projected 
into the {pr,r)-plane after each time interval of length 2tt/VL. The 
figure suggests that dissociation will occur for values of e that exceed 
e ^ 0.002. 

To gain some insight into the absorption of gravitational waves by 
a Newtonian binary, we have also tested the "rate of dissociation", 
defined to be inversely proportional to the number of iterates of the 
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Poincare map required before the osculating conic of the perturbed 
elhpse becomes a hyperbola, by numerical integration of the Hamilton- 
ian system for various values of the frequency Q and phase shift p of 
the incident wave. We assume here that a = (3; moreover, the initial 
elliptical motion is counterclockwise. Although these experiments are 
somewhat difficult to interpret, one fact seems to emerge regarding the 
sensitivity of the rate of dissociation to the polarization of the wave. 
For fixed fl, the maximal dissociation rate is in the vicinity of p = —tt/2 
while the minimal dissociation rate is in the vicinity of p = 7r/2. This 
rate also depends on f2, but in a seemingly unpredictable manner. 



7. Appendix A: C and S in terms of Delaunay elements 

The purpose of this appendix is to express C = cos 29 and S = 
sin 26 in terms of Delaunay elements. It follows from the relation 
9 = V + g that 

C =r^ cos(25f + 2v) = cos 2v cos 2g — sin 2v sin 2g, 
S =r^ sin(2(7 + 2v) = sin 2v cos 2g + cos 2v sin 2g. 

Moreover, 

cos u — e I sin u 

cosf = z smf = V 1 — ■ 



1 — e cos u 1 — e cos u 

these relations follow from (P), and the fact that by definition v u 
as e ^ (cf. P, Vol. 1, p. 100]). Therefore, 

1 

r cos2v =a {-e — 2ecosu + -(2 — e ) cos2m), 
r ^ sin 2t; =a^ Vl — (sin 2u — 2e sin u) . 

There are classical expansions for cos ju and sinjw in Fourier series 
whose coefficients are expressible in terms of the Bessel function J^, of 
order u. Here, this Bessel function is most conveniently defined by 

Ju(x) := — / cos( ut — X sin t) dt. 
2-K Jo 

Following, for example, J. Kovalevsky p. 49], one finds that 

e °° 1 

cosM = - - + ^ -[J;,_i(z/e) - J,,+i(z/e)] cosui, 

OO 

esinw =2 —Jyive] sin vi, 
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and, for j > 1, 

cos ju = ^ — [Ji,_j(z/e) — J^^ji^ue)] cosui, 

sin ju = ^ — \Jy^j{ve) + Jyj^j{ve)\ sin vl. 

Using these expansions, we obtain the Fourier series given in (11) 
where 

A, =^((2 - e2)[J,_2(z^e) - J,+2(^^e)] - 2e[J,_i(z/e) - J.+i(z/e)]), 

J-^J\-e^ {\^Jv-2[vt) + J.+2(z/e)] - 2 J,(i/e)) . (31) 

We note that the functions C and S are analytic and 2-7? periodic in 
the angle variables I and Moreover, partial derivatives with respect 
to the Delaunay elements can be obtained by differentiation of their 
Fourier series term by term. 

To simplify the expressions for the Fourier coefficients computed 
above, we will use the following elementary identities for the Bessel 
functions p. 48]: 

J'M=\\J.-i{x)~J,^r{x)l (32) 

=^-^J^-i{x) - 2J,(x) + J.+2(x)], (33) 

and Bessel's equation 

x^J'lix) + xJ'^ix) + (x^ - z/2) J,(x) = 0. (34) 

The final expressions for Ay and given in (12) are obtained from (31) 
using (32)- ([3^); in fact, the formula for Ay is obtained by standard 
methods using (32), and the formula for By is derived from the origi- 
nal expressions (31) and (33) after noticing that By is proportional to 
Jy{ye) and using Bessel's equation (pi]). 

8. Appendix B: A binary influenced by a distant massive 

THIRD body 

The purpose of this appendix is to explore the possibility of applying 
our results to the three body problem. For a binary system infiuenced 
by a distant massive third body, our "tidal" approach results in a lim- 
iting case of the celebrated problem of three bodies and the question 
is whether the continuation theory of § ^ would be applicable in this 



GRAVITATIONAL IONIZATION 



33 



case. The existence of periodic orbits in the three-body problem has 
been estabhshed in the classical work of Poincare . 



We study the effect of a massive body, metaphorically the Sun, on 
a binary, metaphorically the Earth-Moon system, where the Sun is 
viewed as giving rise to a periodic perturbation of the Earth-Moon 
orbit by tidal forces. To derive the equations of motion that will be 
considered in this appendix, we will model the Earth- Moon-Sun system 
according to the following scenario. The motion of the Sun is neglected 
due to its great mass, its gravitational attraction brings about the 
motion of the Earth-Moon system as a whole on an almost Keplerian 
orbit about the Sun and its tidal influence perturbs the orbit of the 
Moon about the Earth. It is this latter motion that constitutes the 
lunar problem under investigation here. 

To obtain the mathematical model, let us consider the equations of 
motion of mi and m2 — the Earth and Moon in our approximation, 
respectively — as given by (1), with a single perturbing body, namely 
the Sun, with potential 

*(X) = (35) 

It is useful to transform the equations of motion of mi and m2 from 
Xi(t) and X2(t) to the relative coordinates r = Xi— X2 and the center- 
of-mass coordinates Xcm = (^iXi -|- m2X2)/(mi -|- 7712). In terms of 
these new coordinates, the solar gravitational attraction involves the 
relative coordinates r and R = Xcm — Xq. Let us further assume 
that /i = r/i? <^ 1, where r and R denote the magnitudes of the 
corresponding vectors, and consider the expansion of the solar influence 
in these equations in terms of /i. Using the fact that 

where 77 is a constant parameter with \ri\ < 1, the equations of motion 
reduce to 

%^ = -GoMo^ + 0(/.^), (36) 

+ ^ = - K.,r^ + 0{^,% (37) 
where we have introduced the "tidal matrix" Kij such that 

^« = ^(*«-3^). (38) 
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We consider the external body (the Sun) to be so massive (Mq >> 
777.1,777.2) that it is essentially unaffected by the presence of nii and m2 
(the Earth and Moon, respectively). Thus, we can take X0 = and 
therefore the Sun remains fixed at the origin of the inertial coordinate 
system under consideration. Neglecting terms of order r'^ / B? in (36), 
this equation reduces to the Newtonian two-body equation for the rel- 
ative motion of the centcr-of-mass about the Sun. We take this orbit 
to be slightly elliptic (for example, for the Earth-Moon orbit about the 
Sun, the eccentricity ei is approximately 0.017). The resulting expres- 
sion for R can be substituted into (37) to give the equations describing 
the dynamics of the Earth-Moon system in the presence of the Sun. 
We further assume that the relative orbit as well as the center-of-mass 
motion occurs in the equatorial plane of the Sun. This should be a 
reasonable approximation as the Earth-Moon orbital plane makes an 
angle of approximately 5° with the ecliptic (the ecliptic is essentially 
the plane of the Earth's orbit around the Sun) while the ecliptic makes 
an angle of approximately 7° with the equatorial plane of the Sun. It is 
clear that this "tidal" approach to the three-body problem is somewhat 
different from the standard "restricted" approach; in the latter case, 
the mass of the Moon is effectively set equal to zero. 

The Earth-Moon orbit about the Sun has a small eccentricity; there- 
fore, the tidal matrix in (37) will be determined to first order in the 
eccentricity. To this end, let = GoMq/uq (with a© being the 
semimajor axis of the Earth-Moon orbit around the Sun) and note 
that the eccentric anomaly is -u ~ f2t + ei sin Qt, the true anomaly is 
V flt + 2ei sin fit, and i? a; 00(1 — 61 cos fit). Using R} = R cos v and 
B? — Rsinv, the Cartesian components of the tidal matrix are given 

by 

13 3 

Ku ^-fl'^{- + - cos 2 fit - -eiC0sflt{3 - 7 cos2flt)), 

Zi Zi z 

3 

Kx2 = Q^(sin2nt + Ci sin Ot(3 -|- 7cos2ni)), 

z 

13 3 
= - fl^{- - -cos2Qt+ -eiC0sQt(5 - 7cos2Qt)), 

2 2 2 ^3g^ 

i^l3 =i^23 = 0. 

As the tidal matrix is traceless and symmetric, the above equations 
determine all of its elements. 

Using (37) and (39) we can write the associated Hamiltonian for 
this system. Because the motion is taken to be in the equatorial plane, 
polar coordinates are convenient. In these coordinates the Hamiltonian 
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IS 

+ 3 cos 26* [cos 2Vtt - d cos fit(4 - 7 cos 2^^)] 

+ 3 sin 2^ [sin 29.1 + d sin fit(3 + 7 cos 2fit)] }. (40) 

Upon expressing the Hamiltonian equations in terms of intrinsic di- 
mensionless quantities, it becomes clear that the strength of the inter- 
action between the binary and the third body is VL^ /uo'^ <^ 1, but the 
square root of this perturbation parameter also occurs in the harmonic 
terms that render the Hamiltonian (40) explicitly time- dependent. In 
particular, the period of the harmonic terms becomes unbounded as 
the perturbation parameter goes to zero. Therefore, the continuation 
method of § ^ is not directly applicable in this case; in fact, the reso- 



lution of this problem is due to Hill (cf. |^ 1131) • Hill's approach, 
the equation of relative motion (37) is referred to a Cartesian system of 
coordinates r' that rotates with frequency VL with respect to the inertial 
system. Let r* = Sijr'-' , where the nonzero elements of the orthogonal 
matrix S are given by 

5*11 = 5*22 = cos fit, -5*12 = 5*21 = sin fit, 5*33 = !; 

then, the equations of motion in the new system are (r' = r) 

nn^y' n2 / , ^x' 



2Q^ - Q'x' + ^=- {K[,x' + K[,y') 
df^ ' ~" dt " ' r'3 



df^ dt 

^^y' nr^^^' o2 / %' 



+ 2fi— - n'y' + ^ = - {K[,x' + K',,y'), (41) 



where K' = S'^KS, i.e., 

= -2fi2(l + 3eicosfit), fs:;2 = -6fi^eisinfit, K!^2 = --K[-^, 

2 (42) 

and = K!^^ = 0. The system (41) is autonomous for ei = 0. In this 
case, periodic solutions exist as originally demonstrated by Hill and 
Poincare (cf. |jl2|). The continuation of such solutions using ei, < 
ei -C 1, as the expansion parameter can be proved, using the Implicit 
Function Theorem, as originally conceived by Poincare (cf. [|ll| ||^ [|12]) 
for the restricted three-body problem. Of course, the method of § ^ 
is also applicable by following the ideas for isoenergetic reduction as 
discussed in § |. 
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